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Ions accelerated by electric fields (so-called runaway ions) in plasmas may explain 
observations in solar flares and fusion experiments, however limitations of previous 
analytic work have prevented definite conclusions. In this work we describe a numer¬ 
ical solver of the 2D non-relativistic linearized Fokker-Planck equation for ions. It 
solves the initial value problem in velocity space with a spectral-Eulerian discretiza¬ 
tion scheme, allowing arbitrary plasma composition and time-varying electric fields 
and background plasma parameters. The numerical ion distribution function is then 
used to consider the conditions for runaway ion acceleration in solar flares and toka- 
mak plasmas. Typical time scales and electric fields required for ion acceleration are 
determined for various plasma compositions, ion species and temperatures, and the 
potential for excitation of toroidal Alfven eigenmodes during tokamak disruptions is 
considered. 


I. INTRODUCTION 


The phenomenon of particle runaway in a plasma is well known, occurring in both space 
and laboratory plasmas [1 2], It arises because the friction force experienced by a charged 
particle decreases with particle energy, so that the presence of a sufficiently strong induced 


electric field can allow t 
Electron runaway 


;he particle to be accelerated - or runaway - to high energy. 

3| has been extensively studied in the context of magnetic confinement 
fusion in tokamaks, where it can lead to the formation of localized high-energy beams which 
must be carefully controlled {4] . The standard analytic method used to determine the initial 
evolution of the electron distribution function in a fully ionized plasma was introduced by 
Kruskal and Bernstein Idj and later generalized by Connor and Hastie [fjj]. It takes the form 
of an asymptotic expansion of the electron kinetic equation in the electric field strength. Once 
a fast electron population - known as the primary distribution - is established, it can rapidly 
produce further fast electrons through large-angle, or knock-on , collisions [7]. This avalanche 
process leads to the so-called secondary runaway-electron generation, which is often dominant. 
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Ion runaway has long been of interest in the astrophysical community, where it is thought 


Q. 


It has also been 


to contribute to the observed abundance of high energy ions in solar flares 
used to study the behavior of lightning channels [gj and was observed in laboratory plasmas, e.g . 
in the Mega Ampere Spherical Tokamak (MAST) [ 9 J and in the Madison Symmetric Torus 10]. 
The detailed mechanism of ion runaway differs from that of electron runaway. Friction with the 
electrons, which are also drifting in the electric held, acts to cancel a portion of the accelerating 
held. In the ion rest frame, the electrons have net motion anti-parallel to the electric held, and 
a test ion will experience two main forces beyond friction against the background of charged 
particles: acceleration in the electric held, and friction due to the electron drift. These forces 
scale differently with ion charge, and the dominant force is either electron friction - with the 
consequence that the ions are dragged along with the electrons - or acceleration by the electric 
held. In a pure plasma the cancellation of electric held acceleration and electron friction is 
complete, but the presence of impurities, neutrals or effects such as particle trapping in a 
non-uniform confining magnetic held allow a finite effective held to remain 0 Q. 


Furth and Rutherford 


12] generalized earlier treatments by adopting a similar expansion 


procedure to that used to study electron runaway. They determined the steady state ion 
distribution function in conditions typical of operational fusion plasmas. Helander et al. [tj 
then considered the initial value problem resulting from the onset of an accelerating electric 
held, produced for example by a plasma instability. An analytic solution for the initial time 
evolution of the accelerated ion distribution function was developed, but it was noted that its 
application was limited due to the low electric helds required for it to be valid. Both of these 
ion runaway studies considered only the presence of a trace impurity population, consistent 
with typical operating conditions in fusion plasmas. 

Plasmas with significant impurity content are also common, however. Astrophysical plas¬ 


mas often consist of a mixture of dominant species, as well as containing trace elements l], ]13] 


Disruptions [2j||, which are a common cause of electron runaway in tokamaks, are also typically 
associated with an increase in impurity content, either due to deliberate gas injection for mit¬ 
igation purposes or due to plasma-wall interaction. Therefore, in Ref. |14(, the initial value 
formulation of the problem posed in Ref. [ 9 J was extended to account for arbitrary plasma 
composition. The potential for ions accelerated during a disruption to excite low frequency 
plasma instabilities was considered analytically. The results were inconclusive since the asymp¬ 
totic expansion used to develop the analytic solution was not strictly valid for disruption-type 
parameters. The limitations of the analytic solutions available in previous work motivate the 










3 


development of a numerical tool to allow detailed study of the time evolution of an ion runaway 
distribution. 

Here we describe the formulation and implementation of an efficient finite-difference- 
spectral-method tool, CODION (COllisional Distribution of IONs), that solves the two- 
dimensional momentum space ion kinetic equation in a homogeneous plasma. CODION solves 
the ion Fokker-Planck equation as an initial value problem and allows for time-variation of the 
electric field and bulk distribution parameters (temperature, density, charge and mass) of each 
plasma species independently. Due to its speed it is highly suitable for coupling within more 
expensive calculations, e.g. studies of instabilities driven by the fast ions or comprehensive 
modeling of ion acceleration with self-consistent coupling to solvers of Maxwell’s equations. 
Using CODION we obtain illustrative two-dimensional ion velocity space distributions, which 
demonstrate the typical behaviour of runaway ions in a variety of physical scenarios. We show 
that during typical tokamak disruptions ions are unlikely to be accelerated to velocities high 
enough for resonant interaction with toroidal Alfven eigenmodes (TAE). Therefore the experi¬ 
mentally observed TAE activity cannot be explained by the ion runaway mechanism. 

The rest of the paper is organized as follows. In Section [TT| we describe the ion Fokker- 
Planck equation, and in Section Hill we outline its numerical implementation in CODION. In 
Section m we explore the numerical solution, including the effect of various approximations 
to the collision operator. In Section IVl the model is applied to a variety of physical scenarios, 
illustrating typical acceleration time scales in laboratory and space plasmas. Finally we close 
with concluding remarks in Section IVII 

II. RUNAWAY ION DISTRIBUTION 

We consider the problem of ion acceleration by induced electric fields with a component 
parallel to the background magnetic held in a plasma. We restrict ourselves to straight held 
line geometry, and assume a homogeneous background plasma. The time evolution of the 
ion distribution is then given by the Fokker-Planck equation, and particle acceleration will be 
opposed by various friction forces. The effect of friction with neutral particles can be significant 
in various physical situations, for example giving rise to charge exchange losses, which was 
studied in the context of lightning discharges in Ref. [sj. Here we will focus on fully ionized 
quasi-neutral plasmas, in which case the friction is the result of inter-species Coulomb collisions 
only. 

We are interested in the initial value problem where an electric held appears in what 
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was previously an equilibrium state. Therefore we assume that the initial particle distribution 
functions are stationary Maxwellians /m - possibly at different temperatures - and consider 
their distortion from this state by the electric field. We linearize the collision operator about 
this background Maxwellian, and neglect the non-linear contribution to the evolution. This 
restricts the study to situations where only a small fraction of the ion population is accelerated, 
or to the initial stages of ion runaway. A reason for why this is sufficient in the present context 
is that, once a high energy population forms, the runaway ions have the potential to excite 


instabilities 


14], which will have a strong impact on the further evolution of the distribution. 


Note that the non-linear terms of the kinetic equation are sometimes required to account for 
the transfer of energy from the electric field into heating the distribution. These aspects of the 
longer term evolution of the distribution are beyond the scope of the work presented here. 

The linearized collision operator consists of the test and field-particle terms. Momentum 
and energy conservation in ion self-collisions is retained by approximating the field-particle piece 
with restoring terms, as described in Ref. [15(], without the need to include the full field-particle 
operator which is numerically more expensive. Collisions with the other ion species could be 
treated similarly, however this would require the simultaneous evolution of the distribution 
functions of multiple species. Therefore, we consider only the evolution of the ion species with 
the highest runaway rate, so that the other ion species remain approximately stationary and 
only the test-particle piece of the unlike ion collision operator needs to be retained. While it 
is difficult to verify a priori that this condition is satisfied, it can be determined by numerical 
simulation of each ion species individually, assuming the others to remain stationary. Because 
of the sensitivity of acceleration rate to ion charge and mass (as demonstrated in Section |V]) , 
the condition can typically be well satisfied as the acceleration rate of different species is often 
separated by several orders of magnitude. 

The velocity-dependent friction on a test particle resulting from collisions with a 
Maxwellian background species has a peak near the thermal velocity of the background due 
to the form of the Coulomb interaction between charged particles. In the case of electrons, 
the friction force will be a monotonically decreasing function for velocities above the electron 
thermal velocity, al 
synchrotron effects 


owing an electron to run away to large energy (where relativistic p] and 


m We focus on situations where the ion, i, and elec¬ 

tron, e, thermal velocities satisfy vr e Ut*, meaning that their temperatures are sufficiently 
similar that T e /Tj y/m e /mi. Then, if an ion is accelerated away from the bulk, friction 
against the electron population will increase with velocity. This has the consequence that the 
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ion acceleration will be naturally balanced by the electron friction 


electric field is below a threshold value similar to the Dreicer held 


or some v < vr e , if the 


17] for electrons. 


Since an electron reacts to the electric held on a time-scale m e /mj times shorter than the 
ions, we assume the electron population to be in a quasi-steady state on all time-scales of interest 
for ion acceleration. Parallel force balance for the electron distribution then requires that the 
total electron-ion friction cancels the acceleration by the electric held, n e eE\\ = Ylj Rej\\ = 
Zjf dv m e v\\C e j, where the sum is taken over all ion species j in the plasma. Due to the 
small mass ratio, the electron-ion interaction is dominated by pitch-angle scattering, so that 
C e j = C e irijZj /riiZf, and we can solve for the friction between electrons and the ion species of 
interest, R ei y = (n^Zf/Z cf[ )eE y, where the effective charge is = Ej n jZj/n e . 

The electron distribution can be written as f e = fue + fe i> with f^e a Maxwellian 
distribution drifting with the bulk ion velocity V,, and a correction f e \ varying over velocities 
of order Vx e , accounting for the electron drift behavior in the electric held. Then the linearized 
ion-electron collision operator, neglecting terms quadratic in f e i, can be simplihcd 
that momentum conservation in binary collisions requires that R; e = —R e ;, 


18], noting 


C ie {fi,fe} = 


R r 


dfi 


+ C ie {f i ,f Me ( v — Vj)} 


( 1 ) 


mifii dv 

The hrst term, which describes the friction arising from the drifting electron distribution and 
was calculated above, readily combines with that describing acceleration by an electric held, 
giving rise to the effective electric held E* = E\\ — R ei \\/niZie = (1 — Zi/Z eS )E^. Thus, as noted 
in the Introduction, in a pure plasma where Z % = Z e g, net ion acceleration will not occur. Light 
ions with Z, < Z e g can be accelerated in the direction of the electric held. Heavy impurities 


with Z % > Z c ff will be accelerated in the opposite c 
runaway (the latter case was studied by Gurevich 


irection, that is in the direction of electron 


111]). The second term in Eq. (JT|) describes 


the slowing down of the fast ions on the electrons, as well as the slow energy exchange between 
the bulk species. Note that the ion how velocity correction is time dependent and formally small 
in the runaway density. This term will become more significant in the ion-electron momentum 
exchange as the runaway distribution builds up. 

Finally then, the kinetic equation we consider for the evolution of the ion distribution 
function in the presence of an accelerating electric held and arbitrary plasma composition is 
the following: 




+ — E* + 


dt 


rrii 


d i-ed 


' dv 


dt, 




( 2 ) 


where £ = vn/v, and the effect of collisions with the background Maxwellian populations is 
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described by the sum over all particle species s in the plasma, 


\ ' ( • j ; i _ 1 n s Zs f ( x s) — G (x s ) d 

^ %S 1 T ie ^ n e 


]__d_ 

x ? dxi 


2x\ 
„ 2 / 




(i ^ f) f 

2 Tjrx'fG (x s ) fi + XiG (x s ) 

E dxi 




[2u s {v)x^u\\ + v E (v)xlQ] f M i , 


( 3 ) 


where r” 1 = n s e A Zf Z'^ hi A./Anx s = v/vt s = \/m s u 2 /2T s and the usual error function 
4>{x) = (2/^/n) f^dye~ y2 and Chandrasekhar function G(x) = [0(x) — x(p'(x)] /2x 2 appear. 
The dimensionless moments u \i and Q of the distribution function appearing in the momentum 
and energy restoring terms of the self-collision operator are 


3 / d 3 u v s {y)v\\fi 

U\\{fi} = . 2 , 

2 J d i vu s {v)v 2 f M i 


Q{fi} = v Ti 


2 / d 3 vu E (v)v 2 fi 


J d 3 vv E (v)v 4 f M i 


( 4 ) 


with the scattering frequencies given by 


vJv) = 4- 


Xi 


v E (y) =2(4 


G(xi) 4>(xi 


Xi 


xi 


( 5 ) 


The runaway behavior of interest can be demonstrated by considering the simpler dynam¬ 
ics of an ion test particle js, 19] in the presence of the electric held. The ion equation of motion 
takes the form mjdu||/dt = ZieE* + iR est , where the test particle collisional friction is given by 


F 


test / 


V = 


-Z?eE D J2 nsZ * Tt 


rip 


E 


1 + — ) G(x s 
rrii 


( 6 ) 


and E e = n e e 3 In A/Ane^E is the Dreicer held. Thus ion acceleration can occur when E/Ed > 
_F. test jZ\eE d |1 — Zi/Z e fj|. The required electric held values are illustrated in Fig. [H for low 
and high effective charge. The hgures illustrate the non-monotonic ion friction force, with one 
maximum near the thermal velocities of the ion species, and another near the electron thermal 
velocity. Therefore, as hrst described by Furth and Rutherford 12J, for a sufficiently strong 
electric held we may expect ions to be accelerated from their initial velocity to a higher velocity 
at which friction on the electrons dominates, giving rise to a suprathermal population in the 
plasma. We will compare the characteristics of the test-particle behavior, governed by the 
friction force illustrated in Fig. |T] to the numerical solution of Eq. (J2j) in the next section. 

Figure Q] also illustrates that the electric helds needed to accelerate ions are highly sen¬ 
sitive to ion charge and plasma composition, due to their effect on the effective electric held 
E* = (1 — Zi/Z e g)E. Note that the electric helds needed to accelerate ions beyond the electron 
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FIG. 1: Electric field needed to accelerate fully ionized test particles of various ion species in an impure 
deuterium plasma, as a function of velocity, (a) rq/^d = 0.4%, nne/wD = 5%, Z e g = 1.2 and (b) 
uc/u-d = 4%, u-He/uD = 5%, Z e g = 2, and all particle species are taken to be at the same temperature. 
Since Z e g = Zn e in (b), no electric field for which the model is valid will accelerate helium ions. The 
quantities shown are independent of the density and temperature of the plasma. 


thermal velocity are significantly larger than the minimum electric held necessary for accel¬ 
eration. Such strong fields will not be considered here, since the validity of the linearization 
typically breaks down before the ions reach a significant fraction of the electron thermal velocity. 


III. CODION 


In this section, we outline the implementation of Eq. (J2J) in the numerical tool CODION, 
which solves the ion Fokker-Planck equation numerically as an initial value problem to give the 
evolution of the ion distribution function in the presence of an accelerating electric held. It 
uses a continuum-spectral discretization scheme based on that used in CODE 20]. Illustrative 
solutions for a tokamak-like plasma are presented, and a comparison of the obtained distribution 
function is made with the behavior predicted by test-particle equations, demonstrating the 
importance of collisional diffusion for the runaway of ions. In addition we investigate the effect 
of different choices for the self-collision held-particle operator. 

The pitch-angle dependence of the distribution function is represented by a truncated 
Legendre polynomial expansion, while velocity is discretized on a uniform grid v = v n = nAv, 
















n = 0, 1, N — 1: 


fi(v n , 6 , t) = Y2fl(v n , t)Pi (£), 


1=0 


(7) 


where the Legendre polynomials Pi obey the orthogonality relation d£ P;(£)P//(£) = 
(2/ + 1), and 


fl(Vn, t) 


21 + 1 
2 


d zm)fi(v n , e, t). 


'-i 


( 8 ) 


The integral operation (L + 1/2) f_ ± d£ P/-, (£)... is applied to the kinetic equation in Eq. (J2]) for 
each L, producing a linear set of equations for the quantities /l(u, t), using the boundary con¬ 
dition fi max (v, t) = 0 for all v. Well-known recurrence relations for the Legendre polynomials 
are used to obtain analytic expressions for all the terms appearing in the equation. For exam¬ 
ple, the Legendre polynomials are eigenfunctions of the linearized collision operator, while the 
electric field-term will produce a coupling between and fi ±i modes. This procedure exactly 
captures number conservation for any choice of Z max > 1. The velocity derivatives appearing in 
the kinetic equation are represented with five-point stencils 


d A 

dv 

d 2 f, 


1 


N -1 


dv 2 


^ ^ ( &n,m —2 + &&n,m—l 8<5 nim _|-i T 8 n m ^_ 2 ) (Wm), 

m =0 
N—l 

12Au 2 ^ 1 ( ^n,m— 2 T 30(f njm T 1 H$n,m+1 ^n,m+ 2 ) fl(+n 


12Au 

1 


(9) 

( 10 ) 


m =0 


formally introducing an error of order 0( Au 4 ). The integral moments of the ion distribution 
appearing in the self-collision model operator are discretized with a quadrature of the form 
f dv A(v) = w mA(v m ), where the quadrature weights w m are chosen according to Simpson’s 
rule, also with error of order 0( An 4 ). 

For the distribution function to be single-valued and smooth at the origin, we enforce the 
boundary condition f t (0) = 0 for all l > 0 and df 0 /dv\ v=0 = 0. Since we restrict ourselves 
to cases where electron friction will dominate the electric field at sufficiently high velocities, 
the maximum resolved velocity can always be chosen so that only insignificant numbers of 
particles are near the edge of the grid, minimizing the effect of the choice of boundary condition. 
We use the Dirichlet boundary condition fiiv^) = 0 for all l at the maximum velocity. A 
detailed investigation of the convergence properties of the solutions with respect to discretization 
parameters is described in Ref. 21]- 
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The discretization procedure outlined above casts the kinetic equation, Eq. (J2]) , into the 

form 

dF 

—+ MF = 0, (11) 

where M is a sparse matrix and F is a vector representing the discretized distribution function. 
Time integration is performed with the first order implicit scheme 

F(t n+1 ) = [I + A tM[t n+l )}~ 1 F(t n ), (12) 


where any time-dependence of the operator M is due to time-variation of electric fields and 
background plasma parameters. An arbitrary plasma composition is determined by a set of 
input vectors containing particle masses m s , corresponding charge states Z s , charge densities 
p s = n s Z s and temperatures T s . 

Figure [21(a) shows a typical example of the evolution of the ion distribution for a case 
where the electric field is above the minimum required for runaway acceleration. The plasma 
parameters used are characteristic for tokamaks with a hot bulk deuterium plasma at 1 keV and 
fully ionized native carbon impurities. Note that the loop voltage is typically < 1 V in normal 
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FIG. 2: Deuterium distribution function in a hot plasma characterized by T = 1 keV for all particle 
species, n e = 3 • 10 19 m~ 3 and effective charge Z e $ = 2 due to fully ionized carbon impurities with 
nc/nD = 4%. The electric field E = 1.64V/m corresponds to E/Ed = 0.13. (a) Time evolution of 
the | = 1 cut of the distribution and (b) contour plot of the steady state distribution, established after 
20 s. 
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0,0 


tokamak operation, corresponding to E < 0.2 V/m [2|, .9]. A contour plot of the distribution 
in velocity space when steady state is reached is shown in Fig. [2](b). For the parameters used 
here, approximately 5% of the ion population has been accelerated and the linearization used 
to derive Eq. (J2J) is well satisfied. 

The steady state distribution is typically established in 10-20 s at this temperature and 
density, and the time-scale varies with plasma parameters as the collision time defined in con¬ 
nection with Eq. (J3J), r ie oc T^ 2 /n e . For stronger electric fields, the initial evolution of the 
distribution can be followed, but the linearization breaks down before the steady state can be 
reached. The entire ion distribution will eventually run away when E > 0.2 Ed for the Z e g- = 2 


case considered here, and the linearization breaks down within ~ 30 ms with such an electric 
field. 


IV. RESULTS 

Expansions of the collision operators appearing in Eq. (J2J) have been used previously in 
the literature to consider ion runaway analytically. Here we compare the effect of various ap¬ 
proximations to the collision operator on the numerical solution of the distribution function. 
The plasma parameters of Fig. [2] are used as a basis in the comparisons, but the behavior 
illustrated is characteristic of a wide range of parameters. We first consider the effect of ne¬ 
glecting the conservative field-particle terms in the self-collision operator. Figure [3] shows the 
£ = 1 cut of the distribution of the bulk deuterium population for two cases, with effective 
charge Z e s = 1.5 and Z e g = 2, respectively. It can be seen that, using only the test-particle 
operator, the dominant behavior of the fast-ion population as given by the fully conserving case 
is reproduced, with the main difference being the rate at which it builds up. This is expected, 
since the conserving field-particle terms are proportional to / m * oc exp(— xj), and therefore act 
mainly on the thermal bulk of the distribution. 

Figure Ufa) shows a case with Z e g = 1.5. With a lower amount of impurities to which 
momentum and energy can be transferred, the fully conserving linearized collision operator for 
self-collisions will exhibit unphysical behavior before a significant runaway population has time 
to form, which is clearly illustrated by the distortion near v — 0. The reason is that when 
the conserving terms are kept in the kinetic equation, the distribution is heated by the electric 
field, causing the linearized equation to break down after some time. This is also observed 
in the solution obtained using only the energy conserving term, albeit less pronounced. The 
distribution functions obtained with the momentum conserving self-collision operator tends to 
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FIG. 3: Comparison of the £ = 1 cut of deuterium distributions calculated using CODION, retain¬ 
ing all (solid) or individual (dashed, dash-dotted, dotted) conservative terms in the ion self-collision 
operator. Here, T = 1 keV for all species, and ri£)+ = 3 • 10 19 m -3 . The only impurity is fully ionized 
carbon, of density such that the specified effective charge is obtained, and electric fields and times are 
chosen to produce a significant runaway population, (a) Z e g- = 1.5, E = 2.5 V/m at t = 0.9s, and (b) 
Z e g = 2, E = 2 V/m at t = 0.8 s. 


stay regular for longer. Figure |3](b) shows a similar case but with higher impurity content, 
Z c r = 2, for which all operators yield well-behaving results. The main difference between 
the conserving operators and only using the test-particle operator is that the former typically 
lead to a runaway rate which is at least twice as large. An additional consequence of the low 
effective charge was demonstrated in Fig. QJa), where impurity ions were shown to be more 
easily accelerated by an electric field than the bulk species, implying that for low Z e g the 
assumption of stationary impurity ions may be violated. 

In conclusion, the high-energy part of the ion distribution obtained using only the test- 
particle operator is in qualitative agreement with the result obtained with conservative oper¬ 
ators, but the runaway rate is expected to be lower in the test-particle case. A quantitative 
investigation of runaway rates for impurity species is presented in Section |V] 

It is instructive to compare the behavior of the numerical solution to the characteristic 
behavior indicated by the test particle friction given in Eq. ([6]). Noting that the runaway ion 
velocity satisfies Vn v -C VT e , we can expand the contributions to Eq. ((6]) using the known 
low and high-velocity limits of the Chandrasekhar function. The test particle friction in this 

















12 


limit reduces to 


F: 
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nriiVTi 


-'eff 
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3/2 



X, 


(13) 


Tjp 3-v/7r y T e 

where n = JA njZ?rrii/n e mj allows for arbitrary impurity content. Consider first the minimum 
value of the magnitude of the collisional friction force; this will determine the minimum electric 
held which can accelerate a fast test ion. Minimizing Eq. (1T3|) yields 

1/3 


ntest / 


= V Te 


= -2 


3 m e 
2 rrii 
rriiVTi Tj 
Tie. T f , 


(Z e ff T Tl 
3 m 


11/3 


2tt rrii 


(Z eff + n) 


(14) 

(15) 


From this it follows that the minimum, ’’critical”, value E ci of the electric held above which a 
test ion can be accelerated is given by 

1/3 


Zi (Z ( 


eff 


• n) 1 ' 3 


id 


Zi/Z ( 


eff | 


3 m e 
2n rrii 


(16) 


By taking ZieE* + FF si = 0, we can hnd the range of test ion velocities, v c \ < v < v C 2 , for which 


acceleration in a given electric held occurs, as discussed in Refs, [l, [12J. Using the expression 


for the friction given in Eq. (fT3]l results in a third order equation, however simpler approximate 
formulae can be obtained by noting that u cl will fall near to the region dominated by ion 
friction, and v C 2 in the region dominated by electron friction. Retaining only the corresponding 
terms in Eq. (TT3j) yields, for arbitrary impurity content, 


Vd = v Ti 


ZjT e 

~2T~ 
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-l 


Z c ff + n 


|1 — Zi/Z t 


eff 


3 y/r E |1 - Zi/Z, 


eff 


V C 2 ~ V T e- F 7 

These equations generalize the corresponding expressions in Ref. [1] 


(17) 

(18) 


to arbitrary plasma com¬ 
position. Note that these formulae are only valid when E is sufficiently large compared to E ci , 
since at E = E C i ion and electron friction contribute equally. We may expect that in steady- 
state, the position of the high-velocity maximum of the distribution function, denoted v m , is 
close to v c 2 , which scales linearly with E in the approximate form given by Eq. This is 

confirmed numerically and illustrated in Fig. [4j where we show the variation with electric field 
of v m , obtained from steady-state CODION solutions of Eq. ([2]). Also shown are the boundary 
velocities of the acceleration region, resulting from numerical solution of the force balance us¬ 
ing the full test particle friction, Eq. (|6]), as well as their approximate forms Eqs. (UZHUD. The 
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FIG. 4: Electric field dependence of the location v m of the maximum in the runaway tail, obtained 
using CODION, for a fully ionized impure deuterium plasma with nc/n d = 4%, Z e g- = 2.5 and equal 
temperature for all particle species. The boundary velocities and v* 2 of the acceleration region 
obtained numerically from Eq. ([6]) (using the velocity which minimizes the friction when E < E c i), 
and their approximate values v c \ and v C 2 given by Eqs. (I17H18I) . are also shown. The quantities shown 
are independent of electron density and temperature. For electric fields E < 0.8 Ed, no maximum 
forms in the runaway ion distribution. 


values converge when the system is strongly driven by a large E. The linear dependence of v m 
is clearly visible at large E , where it approaches the value given by the test particle approx¬ 
imation. The analytic approximation for Ed, Eq. (fT6j) . is only accurate to ~ 20%, however, 
indicating that collisional diffusion contributes significantly to the evolution at lower electric 
fields. Since the linearization breaks down more rapidly at larger electric fields, it is mainly at 
fields near the threshold for runaway generation that the model can consistently be applied to 
study the long-term evolution of the ion runaway tail, making a full kinetic simulation essen¬ 
tial for capturing the important physics. For the more massive impurities, the features of the 
test-particle model become increasingly accurate since the runaway ion population is further 
separated in velocity space from the thermal bulk. 

It is important to point out that neither the diffusion terms nor the field-particle self¬ 
collisions have been accounted for in the derivations of the above estimates, which are meant 
to give simple expressions that show how the essential quantities scale with the plasma param¬ 
eters, and to provide a useful physical picture for illustrating the ion runaway phenomenon. A 
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complete description will be provided only by numerical solution of the kinetic equation. 


V. APPLICATIONS 

In this section, CODION is applied to calculate runaway ion distributions for typical 
solar flare and fusion plasmas. The rate at which a fast ion population forms due to the 
runaway mechanism is determined, and it is investigated whether the difference in acceleration 
rate between different ion species can explain the enhanced abundance of heavy ions in the 
solar wind. We also consider the possibility of Alfvenic instabilities being driven by runaway 
ions during tokamak disruptions. 

Throughout this section, time-scales are chosen so that significant fast ion populations 
have time to form, which typically takes a few hundred ion-electron collision times. We define 
the runaway density, n r = (1 /rti) f v>v * d 3 v fi, as the fraction of ions with velocity larger than 
the low-velocity numerical solution of ZieE* F^ est = 0, denoted Ki, which if E < E cl is taken 
as the velocity u m ; n minimizing the friction i7 test . 


Ion acceleration in solar flares 


Solar flares are thought to be initiated by reconnection in the corona |22j, but the origin of 
the observed fast ion populations in flares is still not completely understood |23j. Both stochastic 
acceleration by waves and the direct acceleration of the particles by the electric held have been 
considered, and it appears likely that a combination of the two can be at work 


25]. 


The effective accelerating held experienced by a given species varies with its charge and 
the effective charge of the plasma, as discussed in Section [TT1 This can give rise to preferential 
acceleration of heavier elements under certain circumstances, and this effect was considered 
in II where elates of the runaway rate were gl ve„ baaed on the approximate formula 
of [11]. With CODION we can determine the time evolution of the ion distribution function 
numerically, and evaluate the dependence of the acceleration rate on various ion parameters. 

The composition of the solar plasma, particularly the metallic elements, has been studied 
extensively in recent years, however much uncertainty remains. We will choose parameters 
consistent with the choices made by Holman |lj]. We use the plasma temperature T = 700 eV 
for all particle species, and hydrogen density rift = 3 • 10 1 ' m -3 . Elements with atomic number 
Z < 6 can readily be assumed to be fully ionized at this temperature. The plasma composition is 
based on the ion abundance recommended by Schmeltz et al. [13]. We use a helium population 
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FIG. 5: Distribution functions of hydrogen, helium and carbon after 30 s of acceleration in a solar 
flare-like plasma, with E = 50 mV/m. All distribution functions are normalized to their respective 
densities. The temperature is T = 700eV for all species, nn = 3 • 10 1 ' m" 3 , bhc/bh = 6%, and all 
other impurities are represented by a carbon population with bc/bh = 0.1%, yielding Z e g- = 1.13. 
The runaway densities are n r H ~ 0, n r , He = 3.7 • 10 4 and n ri c = 0.18. 

of density nue/nu = 6%, and represent all heavier impurities by a carbon population of density 
nc/n H = 0.1%, corresponding to an effective charge Z e s = 1.13. Electric field strengths in solar 
flares are not well constrained by experimental observation, and we will investigate the rate of 
acceleration at a range of values. The Dreicer field is Ec = 224 mV/m for this set of plasma 
parameters. 

The critical electric fields E C i for the ion species in such a plasma are E c ^ = 154mV/m 
for hydrogen, E c ^ e = 40 mV/m for helium and E C) q = 20 mV/m for carbon. Note that the 
acceleration rate depends strongly not only on E/E ci , but also on v c i/vn- Figure [5] shows the 
v± = 0 cut of the distribution functions of hydrogen ( 1 H), helium ( 4 He) and carbon ( 12 C) after 
30 s of acceleration from initial Maxwellians, with the plasma parameters specified above and 
an electric field E = 50 mV/m. This is significantly below the hydrogen critical field, and no 
hydrogen runaway population forms. Runaway ion populations of both helium and carbon do 
form however, with runaway densities n r ^ e = 0.037% and n r) c = 18%, respectively. 

Positive values of ny represent the direction of the electric field. Therefore, Fig. [5] also 
illustrates how heavier ions (charge Z > Z eS ) are accelerated in the direction opposite to the 
electric field, dragged by electron friction. The corresponding 2D carbon distribution function 
is shown in Fig. [6l displaying a strong directional anisotropy (compare to Fig. E)(b)). This can 
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FIG. 6: Distribution function fc of carbon, accelerated from an initial Maxwellian for 30 s by an electric 
field E = 50mV/m, in a fully ionized solar flare-type plasma with nn = 3 • 10 17 m -3 , nne/^H = 6% 
and nc/un = 0.1% , and temperature T = 700 eV for all species. 

be understood by the observation that the accumulation velocity v C 2 is located at a higher value 
of Xi = v/vTi for heavier ions. Since pitch angle scattering of the energetic heavy ions scales 
with velocity like Z e s/xf, the mechanism will be less effective than for light bulk ion species in 
increasing the perpendicular energy of the distribution. 

We will now investigate how the rate of acceleration varies with electric field strength for 
different ion species. To determine the runaway density of heavier ions, we have introduced 
trace amounts of each ion species with charge between 2 (helium) and 18 (argon), assumed to 
be fully ionized. In practice, ions of charge Z* > 8 will typically not be fully ionized at the 
temperature considered due to their high ionization energy, meaning that the results shown 
here will overestimate the acceleration rate of the heavier ions. The ion masses have been set 
to that of the most common isotope, i.e. 'Li, 9 Be, 20 Ne etc. Both 3 He and 4 He are shown, with 
3 He showing a significantly enhanced runaway rate as compared to that of 4 He, for all values 
of E. 

Figure Ha) shows how the runaway density n r depends on electric field after 1 s of ac¬ 
celeration for various ion species present in the solar flare plasma. The figure illustrates how 
the runaway rate is sensitive to ion parameters; at low electric fields the heavier ions tend to 
be accelerated at a lower rate than light ions, while at higher fields they are the most readily 
accelerated. Note that above the critical electric field, E cl ps 40 mV/m for helium, the runaway 
rate increases significantly faster than for E < E C i. The runaway rate of 3 He is seen to be 
orders of magnitude higher than the runaway rate of 4 He for all electric fields considered. 

Finally we illustrate the dependence of the acceleration rate on ion charge and mass. 
Figure Hb) shows the runaway density n r after 1 s of acceleration as a function of ion charge Z % 
for various electric fields. Ions of charge between Z* = 4 ( 9 Be) and Zj = 8 ( 16 0) are seen to be 
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FIG. 7: Solar flare-like plasma with all particle species at the same temperature T = 700 eV. The 
composition is mainly hydrogen, nn = 3 • 10 1 ' m~ 3 , and helium-4, nne/^H = 6%, with carbon of 
density nc/nn = 0.1% representing all other impurities, each appearing in trace amounts. Ions are 
accelerated for Is from initial Maxwellians. Due to the low effective charge Z e s = 1.13, the hydrogen 
background remains stationary ( E* is small), and heavier elements are assumed to be present only in 
trace amounts, making self-collisions negligible, a) Runaway density n r as a function of electric field 
strength E for various ion species, b) Runaway density as a function of ion charge Z % . For each atomic 
number the most common isotope is used, except for helium (Zj = 2) where both 3 He and ^He are 
shown, with the 3 He runaway density indicated by dots. 


preferentially accelerated over lighter or heavier elements for low electric fields. For Z, > 8, the 
trend depends on electric held. For low electric fields, the runaway rate decreases with charge, 
while for larger electric fields (in this case above approximately 50 mV/m), heavier ion species 
may be accelerated at a higher rate than the lighter species. Further studies are required to 
determine the full effects of variation in background composition, temperature, density, and 
charge states on the relative rates of acceleration between different ion species. The results 
presented here demonstrate the utility of CODION for the problem. 

As previously noted, acceleration by quasi-static electric fields is not the only mechanism 
for ion acceleration in a hare plasma. Interaction with Alfven waves can accelerate ions which 
have velocities above the Alfven velocity. As this is usually well above the thermal ion velocity, 
an initial acceleration by electric helds may be required before the process becomes significant 1 


26]. CODION provides the means for more accurate modeling of the effects of such interactions. 
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B. Tokamak disruptions 


During tokamak disruptions the plasma temperature drops from the typical operating 
regime of around several keV to a few eV in a couple of milliseconds. A large electric field 
is initially induced parallel to the magnetic field to maintain the plasma current of several 
MA, potentially leading to the formation of a beam of energetic electrons through the runaway 
mechanism. The potential for damage by such a focused high-energy beam on contact with the 
vessel wall is large, and runaway generation must as far as possible be suppressed. To study 
the physics and mitigation of runaway electrons, disruptions can be induced by the injection 
of large quantities of noble gas, often in amounts comparable to the initial plasma inventory or 

larger [2J]. 

The large induced electric field will usually decay rapidly on a timescale of a few ms in 
response to the formation of a narrow runaway electron (RE) beam. With runaway electrons 
reaching energies of order tens of MeV, they can carry a significant fraction of the pre-disruption 
plasma current and can drive high frequency electromagnetic instabilities through resonant 


interactions 


27 


30]. Recently, low-frequency magnetic fluctuations in the range / 


60 


260 kHz, have been observed in the TEXTOR tokamak during induced-disruption studies with 
argon massive gas injection (MGI). These fluctuations take the form of either a strong signal at 
a distinct frequency |31|, or accompanied by broadband activity [32]. The fluctuations appear 
to limit the RE beam formation in these cases, as the magnetic perturbations may scatter 
the runaway electrons and provide passive mitigation. Aside from the potential consequences 
for mitigation, observed instabilities offer a non-intrusive diagnostic for both bulk plasma and 
fast-particle properties, through the extensively applied technique of MHD spectroscopy [331 ]. 


Fast ions resulting, for example, from certain heating schemes are well 


drive low frequency Alfvenic instabilities in typical operational scenarios 3J, [35|. Runaway ions 


mown to resonantly 


may thus also provide a potential drive for the fluctuations observed. Interestingly, toroidal 
Alfven eigenmodes (TAEs) 0 can have frequencies and mode numbers in the same range as 
the post-disruption magnetic activity. Therefore, the excitation of TAEs by runaway ions was 
recently considered in Ref. [dJ] using an analytical approximation for the runaway distribution 
function. As the validity of the approximate distribution was limited, definite conclusions could 
not be drawn. With CODION we can extend the study using the numerically calculated ion 
distribution. The cold post-disruption plasma is highly collisional, motivating the use of a 
homogeneous background plasma and the neglect of magnetic trapping when evaluating the 
effective electric field. 
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The TAE perturbation is typically dominated by two neighboring toroidally coupled har¬ 
monics at large aspect ratio, with poloidal mode numbers m and m +1. The mode is localized 
about the minor radius r = ro at which the magnetic safety factor is go = (2m + l)/2n, where 
Ro is the radius of the magnetic axis and n is the toroidal mode number. The TAE frequency 
is uj = va/( 2q 0 Ro), where va = B / ^Jpop m is the Alfven speed and p m the mass density. The 
two component harmonics allow resonant interaction with particles whose parallel velocity W|| 
satisfies |u||| ~ va /3 or |u||| ~ va■ It was argued in Ref. Q that as the runaway ions accelerate, 
the inverted region of their energy distribution, df/dS > 0 where £ = mp/ 1 j2 is the particle 
energy, can reach the lower Alfven resonance, U|| = va /3 and may drive the TAE. If the radial 
runaway ion profile peaks on axis, the spatial gradient df /dr will give an additional positive 
contribution to the growth rate. Taking parameters characteristic of argon MGI-induced dis¬ 
ruptions, no = 3- 10 19 m -3 , n\ T = O.Itt-d, Zat = 2, and anticipating a native background carbon 
impurity with rtc = 0.08no and Z c = 2 so that Z eS = 1.26, a background ion temperature of 
Ti = 10 eV, toroidal magnetic field B = 2T and major radius Ro = 1.75 m means that the 
resonance condition vu = va /3 requires deuterium ions with velocities v ~ 35utd- 

The electric field required to accelerate bulk ions to the resonant velocity at these low 
temperatures is substantial, varying in response to changes in Z e g but is typically > 0.3 Ed ~ 
100 V/m. Such field strengths are unlikely to occur during a disruption, and they would be 
short-lived if they did 371 ]. Therefore, we conclude that whilst ion runaway may be of interest in 
hot fusion plasmas, runaway ions are unlikely to provide the drive for the observed fluctuations 
during disruptions. 

To quantify the electric field needed for significant ion runaway, we show in Fig. [8] how 
the deuterium distribution evaluated after 2 ms of acceleration from an initial Maxwellian - a 
typical time scale for the induced electric field - varies with (a constant) electric field. The 
parameters are nn = 3 • 10 19 m~ 3 , T = 10 eV and the same plasma composition as before with 
Z e ff = 1.26. It is seen that for electric fields below ~ 200 V/m, no runaway tail tends to form, 
and even with E = 260 V/m the fast ions are far from the resonant velocity near 35 uti> The 
behavior is sensitive to which temperature is chosen for the plasma: increased temperature 
decreases the electric fields needed to accelerate ions, but makes the acceleration timescale 
longer. 

Note that a higher amount of assimilated argon, or a weaker magnetic field, would lead 
to a lower Alfven velocity and TAE frequency, and therefore the runaway ions would reach the 
resonance condition more easily. The higher electron density leads to an increased collision 
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FIG. 8: Deuterium distribution function resulting from various electric fields after t = 2 ms of accelera¬ 
tion from an initial Maxwellian in a TEXTOR MGI post-disruption-like plasma, with no = 3T0 19 m -3 , 
nc = 0.08no, n-Ar = 0.1 rip and T = 10eV for all species, with Zq = Zj± T = 2. The critical held for 
deuterium is E c d = 295 V/m, and the Dreicer held is Ed = 962 V/m. 


frequency, however, so that the runaway ion distribution requires a longer time to form. The 
level of argon absorption during mitigation varies between machines |2[ and the absorption 
appears to decrease for bigger machines. Therefore, although we may expect a higher electric 
field in large tokamaks such as ITER, the higher bulk plasma density and likely lower absorption 
suggests no ion runaway will occur. Note that as the ITER magnetic field is so high, even if 
the assimilated argon was equal to the initial deuterium inventory, the ions would need to reach 
~ 40 Vxi (at 10 eV) to reach the typical va /3 resonance, assuming B = 5T, R 0 = 6 m and 
no = 1 • 10 20 m -3 . 

Finally we note that the model presented in this paper assumes the initial ion distribution 
to be a stationary Maxwellian. Fast ion populations present due to heating schemes in use before 
the disruption may not be completely expelled, and it is not certain that the initial distribution 
is accurately described by a Maxwellian. The dynamics of ion acceleration starting from these 
non-Maxwellian distributions may yield a fast-ion population more readily. Furthermore, since 
plasma conditions change on a time scale of a few collision times during the sudden cooling 
associated with instabilities or disruptions, so-called “hot-tail” runaway generation may occur. 
This has been shown to be an important effect for runaway of electrons, where a seed of 









21 


runaway electrons are provided by fast electrons present before cooling 


. These fast electrons 


are cooled at a slower rate than the low-energy electrons, and may find themselves in the 
runaway region when the plasma has reached its final temperature. However, using CODION 
to investigate the effect of hot-tail ions, it has been concluded that the effect is small for realistic 
fields. The reason is that, if the electric field is not high enough for the ions to overcome the 
friction and become runaways in the first place, all hot-tail ions will necessarily also be slowed 
down. A low-velocity inverted ion population may however form as a result of the cooling 
process even for such low electric fields, with its peak near the velocity which minimizes the 
collisional friction at the final temperature (typically around 6-8 thermal velocities). This is 
still significantly lower than that needed for resonant interaction with Alfven waves. 

One refinement to the model would be the inclusion of “knock-ons”, i.e. large-angle col¬ 
lisions, which have been neglected in the Fokker-Planck equation. It is well known that single 
collisions can change the momenta of the interacting particles significantly, and a runaway ion 
interacting with a bulk ion could cause both to end up in the runaway region. In a situation 
where the electric field is low enough that runaway ions are produced at a low rate through the 
standard acceleration mechanism, knock-on collisions could possibly contribute significantly to 
the runaway generation rate. This has been demonstrated to be the case for electron runaway, 
where this effect drastically affects the rate at which runaways are produced. A simplified 
runaway ion knock-on operator could potentially be constructed from the Boltzmann collision 
integral under the assumption that fast ions accumulate near v C 2 and collide mainly with the 
bulk distribution, since the fast ion distribution is assumed to be a small perturbation in our 
linearized model. However, there are differences between ion runaway and electron runaway 
that suggest that knock-on runaway generation is a less significant effect for ions than for elec¬ 
trons. Since our linearized model restricts the study to cases where E ~ E ci (which is also 
the regime where knock-on generation would be expected to be significant), the accumulation 
velocity near v C 2 will not be significantly larger than the runaway velocity v c \. Because of 
this, collision events where both particles end up in the runaway region will be less frequent. 
This is in contrast to electron runaway, where the electrons have unbounded energy (neglecting 
radiation effects). 

There are applications for knock-on operators other than avalanche generation. It has 
been suggested that fast ion populations due to other sources - for example hot alpha particles 
created in fusion reactions or ions heated by external sources such as neutral beam injection 
(NBI) or radio frequency (RF)-heating - could accelerate bulk impurity ions, which could in 



22 


turn be used for diagnostics [39|, [40]. The suggested collision operator could be implemented in 


CODION, and the time-evolution of bulk impurities solved for in the vicinity of an assumed or 
numerically obtained background of fast ions, however this is outside the scope of the present 
paper. 


VI. CONCLUSIONS 


Electron runaway resulting from the occurrence of a strong electric field in a plasma has 
been the subject of extensive study, and numerical tools exist to simulate the electron dynamics. 
The analytic description of the associated ion acceleration was developed at the same time, but 
its application has been much more limited and is restricted by the various approximations to 
the collision operator which were required. 


We have developed an efficient open-source numerical tool, CODION [4l|, which solves 
the ion Fokker-Planck equation as an initial value problem in a fully ionized plasma. A uniform 
background magnetic field is assumed, along with initially stationary Maxwellian distributions, 
however arbitrary impurity densities and temperatures may be specified. A model operator for 
ion self-collisions based on that used in the gyrokinetic code GS2 |42[] has been employed, satisfy¬ 
ing momentum and energy conservation, non-negative entropy production and self-adjointness. 
A simplified analytical model based on the large mass ratio is used for ion-electron collisions, 
allowing a description of ion-electron friction caused by the perturbation of the electron dis¬ 
tribution due to the electric field. However, we wish to note that our model will break down 
for strong electric fields, as the electrons - which are assumed to be in force balance with the 
electric field and ion friction - will be rapidly accelerated by electric fields approaching the 
Dreicer field. A full description of such scenarios would require the simultaneous evolution of 
the ion and electron distributions, for example by coupling the CODION and CODE 0 
codes. 

The effect of various approximations to the collision operator commonly used in the liter¬ 
ature has been studied numerically. It has been demonstrated that the addition of momentum 
and energy conservation in self-collisions mainly acts to increase the rate at which the fast ion 
population builds up, while the qualitative behavior is largely unaffected. For strong electric 
fields, the test particle description is seen to reproduce the characteristic velocity achieved by 
the runaway population well, and the predicted critical electric field for ion acceleration is 
accurate to ~ 20%. Using the test particle approximation, we derived concise analytic expres¬ 
sions for the critical electric field for ion runaway, Eq. (fTfilh and the typical runaway energy, 
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Eqs. (EHUD, and tested these expressions against direct numerical simulation with CODION. 

The output of CODION is the evolution of the 2D velocity space ion distribution. The 
utility of this has been demonstrated for calculating acceleration rates of ions in solar flare 
plasmas. The rate at which ions are accelerated has been evaluated for a range of ion masses 
and charges for a solar flare scenario based on that considered by Holman [ij , and an exponential 
dependence of acceleration rate with charge for Z > 8 has been illustrated for this scenario. 

Low-frequency instabilities, in the range characteristic of TAE modes, have been observed 
in post-disruption tokamak plasmas, where the disruption was induced by massive gas injection. 
Using CODION we have considered the potential for ions accelerated in the disruption-induced 
electric field to drive such modes resonantly. The post-disruption discharge parameters are 
not well constrained experimentally and simulations of a range of values indicate that ion 
acceleration is possible. The typical maximum ion velocity achieved is too low for resonant 
interaction to occur, however, and the rate of runaway generation is too slow for a significant 
runaway density to be reached in the short-lived electric fields of a typical disruption. 
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